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Abstract 

An observer-based Hamiltonian identification algorithm for quantum 
systems is proposed. For the 2-level case an exponential convergence re- 
sult based on averaging arguments and some relevant transformations is 
provided. The convergence for multi- level cases is discussed using some 
heuristic arguments and the relevance of the method is tested via simula- 
tions. Finally, the robustness issue with respect to non-negligible uncer- 
tainties and experimental noises is also addressed on simulations. 
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M , 

c3 ' 1 Introduction 

g. 

, The ability of coherent light to manipulate molecular systems at the quantum 

scale has been demonstrated both theoretically and experimentally [TJ [121 1301 121 
!3|51llS[TgiliniIISir3IIll23llIll2D]- Many of the procedures, considered in 
j_i ■ this aim, are based on the possibility to perform a large number of experiments 

in a very small time frame. Thus, the output provided by these experiments 
can be used to correct the process and to identify more satisfactory control 
fields HH23H7]. 

The ability to rapidly generate a large amount of quantum dynamics data 
may also be used to extract more information about the possibly unknown 
parameters of the quantum system itself. For each test field (i.e., control), there 
is the possibility of performing many observations for deducing information 
about the system, and this process can often be carried out at a much faster 
rate than the associated numerical simulations of the dynamics. Moreover, the 
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recent advances in laser technology provide the means for generating a very 
large class of test fields for such experiments. 

The rapidly developing theory of quantum parameter estimation has been 
investigated through different approaches. The maximum- likelihood methods 
and the subsequent experiment design techniques provide a first class of results 
in this area [TH] [35] EH [HI E3] ■ The optimal identification techniques via least- 
square criteria's [6] [5] [15] and the map inversion techniques ]57J are some other 
techniques explored in this area. Finally Kalman filtering techniques pj [2"8] 
have been applied to some atomic magnetometery problems. 

However, in general, developing effective identification algorithms is of a 
great interest in this domain. The main concerns in quantum parameter esti- 
mation theory are the presence of local minima's for the optimization problems, 
sensitivity with respect to the experimental uncertainties and noises and finally 
the heavy cost of computations in formerly developed algorithms. 

Before going through the identification and the experiment design problems, 
we need to ensure the identifiability of the system. This issue has been addressed 
in a recent work [15j where sufficient assumptions applying the uniqueness of 
the inversion result are provided in two relevant settings. A brief review of 
an identifiability result needed for the purpose of this paper is given in the 
Appendix [A] The semi-constructive proof in [T5] suggests that a well-chosen 
control laser field, coupling all the eigenstates of the free Hamiltonian, would 
be sufficient to identify the unknown parameters. 

In [10], a state observer for a known quantum system is proposed. This 
observer is then used as a basis for the quantum parameter estimation applying 
an iterative search algorithm. The provided optimization algorithms typically 
converge toward local minima's. 

Here, in the same direction, we provide an observer-based parameter esti- 
mation algorithm based on techniques derived from adaptive control theory. In 
this aim, we will integrate online a generalized observer including the estimators 
for both the unknown state and the unknown parameters of the system. 

In the next section, we will present the suggested algorithm on a simple 
2-level system where the unknown parameter to identify is reduced to a real 
constant 9 multiplying the dipole moment. After presenting the system and 
the estimator, we check the performance of the method by a first simulation 
(Subsection l2.2|l . Subsection 12 . 31 has for goal to explain the special choice of the 
estimator ([3]) and ([4]). Finally in Subsection 12. 4[ we provide a detailed proof of 
a convergence result. 

In Section[3] we extend the observer of Sectionrjjto the general case of an TV- 
level system. The convergence of the estimator is discussed using some formal 
arguments. The efficiency of the technique is then checked on two test cases of 
3 and 4 dimensions. 

Finally, in Section 2] we address the robustness issue with respect to the 
measurement and the control noises and uncertainties. The effect of different 
kinds of uncertainties and noises on the identification result has been checked 
out on simulations for the 2-level case. Similar simulations for the multi-level 
situations give rise to the same kind of robustness. 
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2 The 2-level case 



2.1 The system and its estimator 

In this section, we consider the 2-level system, 

*!(:0-<* + ■*>"•>(:;)■ -(:;)** 

y(t) = (P*(t),*(t)) , (1) 

where y is the system's output, being the population of the first eigenstate of 
H. Here > the dipole moment parameter is supposed to be unknown. The 
goal of this section is to identify this unknown parameter. 

In the density matrix language, this same system can be written as, 

j t p(t) = -i [H + u(t) 6 M , p{t)] , y(t) = Tr [Pp{t)] , (2) 

where p{t) = g C 2x2 is the projection matrix on the wavefunction ^(t). 

We consider the laser field u(t) to be in the resonant regime with respect to 
the natural frequency of the system (being the difference between the eigenvalues 
of the Hamiltonian H, —uj/2 and w/2): 

u(t) = A(t) cos(w£), 

where is a slowly variable modulation of the amplitude. Here, for simplicity 
sakes, we consider a constant amplitude A(t) = A > 0. 

In this paper, we propose the following observer-based parameter estimator, 

j t P = -i[H + u(t)/*, p] + T(y(t) - y(t)) (Pp + pP - 2Tr [Pp] p) , (3) 
V = Tr [Pp] , 

j f § = -i 7 U (t)Tr [P[ M , p]] - y(t)), (4) 

where 7 and T are positive real design parameters. As we will see in Subsec- 
tion [231 in theory, we need the following validity domain in order to ensure the 
convergence toward the true parameter 9: 

T<A6, A6^u, 7 < 9. 

However, the simulations of the next subsection show that these restrictions can 
be much relaxed in practice. In particular, one can use instead of |3]) and (j4|), 
the average estimator of Remark [IJ where the high Bohr frequency is removed 
and that does not require the precise knowledge of the transition frequency. 
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2.2 Simulations 



Before explaining the choice of this identification algorithm and going through 
the technicalities of a convergence proof, let us check the efficiency of the algo- 
rithm in a simulation. Here, we assume the unknown parameter 9, the frequency 
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Figure 1: The first figure shows the convergence of the parameter estimator 9 
towards its true value 1; the second shows the convergence of the error term 
e(t) toward 0. 

u>, the constant T and the laser amplitude A to be all equal to 1. The constant 7 
in (U) is chosen to be 7 = 0.1. The real system p and the estimator p are respec- 
tively initialized at p(0) = * *J and ,5(0) = * *S where <f = (l/V%, l/\/2) T 
and ^0 = (l/v5) 2/-\/5) T . Finally the parameter estimator 9 is initialized at 
9(0) = 1.5. 

The simulations of Figure Q] show the result for the algorithm presented 
above. The simulation time T = 507r represents 25 times the natural period, 
at which the system without control oscillates. As it can be seen the above 
algorithm ensures the convergence of the parameter estimator 9(t) toward the 
unknown parameter 9 = 1. 

2.3 Estimator design 

In this subsection, we will explain the particular choice of the observer-based 
estimator © and 

Starting from the physical system |T]), whenever the parameter 9 is known, 
an intuitive observer (similar to the one proposed in |10p can be given as follows: 

j f V = -i(H + 9 u(i)/i)f + T(y(t) - y(t))P*, (5) 

where y(t) = (P^(t), ^(t)) and T is a positive constant. 
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The non-conservative wave- functional equation ([5]) (||^(t)|| does not remain 
constant equal to 1) can be written in the density matrix formulation as follows: 

j t P = -i[H + 9 u(t)p, p] + T(y(t) - y(t)) (Pp + pP) , (6) 

where p = ty^St' and y = Tr [Pp]. 

Note that unlike the conservative Schrodinger equation @, the observer 
equation Q does not conserve the trace of the density matrix p. However, as 
we will see in Subsection l2.4l enforcing the observer to keep the same geometrical 
structure as the main system simplifies considerably the convergence proof. In 
this aim we define a normalized density matrix given by p — p/Tr[p\. This 
normalized observer state verifies a conservative equation of the following form 

j t P = -i[H + 9 u(t)p, p] + T(y(t) - y(t)) (Pp + pP - 2Tr [Pp] p) , (7) 

where y = Tr [Pp] ■ 

Whenever, the parameter 6 is unknown, one might consider it as a part of 
the system's state to be estimated. In this aim, we replace the parameter 9 
in by to be estimated (so that we obtain ||3J)). The question becomes now 
to provide an evolution equation for 8 ensuring the convergence toward the true 
parameter 9. 

Consider the Lyapunov function 

v(t,p,6) = ±( y (t)-mr + ±(6-ey, (8) 

where 7 is a small enough positive constant. Deriving with respect to time, we 
have: 

jV = (i§(t) u(t)Tr [P[p,p(t)]] - i9 u{t)Tt p(t)]]) e(t) 

+ -9{0(t) -9)- 2r|e(i)| 2 Tr [Pp] (1 - Tr [Pp]), 

1 

where e(t) = y(t) — y(t). Let us take 

j t § = -t 7 u(t)Tr[P[p,p]]e(t). 

Then we have: 

= {i9 u(t) Tr [P[p, p - p]]) - 2r|e(t)| 2 Tr [Pp] (1 - Tr [Pp]), (9) 

While the last term in © is always negative, the situation has no reason to 
be the same for the first one. We will see however that under certain circum- 
stances this first term can be neglected and the Lyapunov function represents a 
decreasing behavior in average. 

The Figure [2] illustrates the evolution of this Lyapunov function for the 
simulations of the last subsection. 
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Figure 2: The positive definite function V defined in © with 7 = .1; the 
function decreases in average and ends up converging toward 0. 



2.4 Convergence analysis 

One has the following convergence result: 
Theorem 1. Consider the 2-levels system describes by 

d 



dt 



p(t) = -i [H + u(t)9p,p{t)\ 
y(t) = Tr[Pp(t)} 



(10) 



where p(t) is the density matrix, H = ^o~ z , p = o~ x , P = (1 + a z )/2. The goal 
is to estimate 9 > and p via the measure y and the knowledge of to. Assume 
that the estimations 8 and p obey the following dynamics (a nonlinear filter of 
V): 



-p = -i[H + 6 «(*)/*, p] + T(y(t) - Tr[Pp}) (Pp + pP - 2Tr[Pp] p) 
j t § = -i 7 u(t)Tr[P[p, P }} (y(t) - Tr[Pp\). 



(11) 



where V and 7 are two positive gains. Assume that u(t) = Acos(ojt) where A is 
a constant amplitude. Assume that the p(0) does not belong to { "f* > + 2* } • 
Then exist e > and r\ > 0, such that, for any design parameters (A,T,"f) 
satisfying 

< T < eA9, < AO < eus, < 7 < e9 (12) 

and any initial conditions (/3o,#o) satisfying (po corresponds to a pure state, po 
is symmetric positive matrix of trace one and Pq = po): 

Tr[(p ~p(0)) 2 ]<r h \6 o -0\<n, 

the estimates (p, 9) converge: 

lim (p(t) - p(t)) = 0, lim 8(t) = 9. 
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Moreover, this convergence is exponential (locally): it is robust to small modeling 
and measurements errors. 

Remark 1. The fact that p(0) must be different from (1 ± a x )/2 is not a se- 
vere limitation in practice since an alternative to ([lip will be the first averaged 
system (| 13[) . It provides a more realistic estimator since it does not depends on 
lo the Bohr transition frequency that is large in general: 



dt AO - 
dt 4 



o v t [y{t)-Tr 



Tr 



O-zS, 



2 7V 



O z £, 







More-over if the laser frequency lo t does not perfectly match the resonance con- 
dition u{t) = A cos(w r (t)), the averaged estimator reads 



d 



dt 
d ~ 



.LO-U T . £, 

4 — o — [<?*,£] 



.AO, 



lA 
. - ——— lr 
V dt 4 



2 + i (y{t) - Tr [a z i] ) (cr,| + £<7, - 2 Tr [<r„£] f 

When \u> — w r \ <C ^40, i/ie second averaged system is then identical to l|15p and 
i/iws convergence is ensured. Thus a precise knowledge of lo is not necessary for 
estimating 0. 

Remark 2. The simulations show that the un-normalized observer equation ([5]) 
( and the parameter estimator found by adapting this equation to the Lyapunov 
function jSJ) ) would be sufficient to ensure the convergence. However, as we 
will see in the proof below, the trace conservation for the normalized equation 
simplifies considerably the analysis since p remains on the Bloch sphere. 

The proof is based on two successive averaging: the first one relies on the 
resonant control u(t), removes the laser frequency lo and yields to (|13|) : the 
second one eliminates the Rabi frequency AO (AO <C lo) and provides (fT5|) . This 
second averaged system is proved to be locally exponentially convergent. This 
implies the exponential convergence of the first averaged system (fT3|) and of the 
original system 



Proof. Notice first that ifTTj) reads 



dt 
d_ 
dt 



P 



^a z + u(t)(T x , p + ^Tr [a z (p - p)} (a z p + pa z 



2Tr [a z p}p) 



iu(t) 



Tr [a y p\ Tr [a z (p - p)} 



2 " 

Since Tr [p(0)] = 1 and Tr [/5 2 (0)] = 1, we have Tr [p(t)} = 1 and Tr [p 2 (t)] = 1. 



Thus p remains a positive matrix associated to a pure state 



i> ) S C 2 of length 



one: p 



p = 



1 + Tr [a x p] o~ x + Tr [a y p] a y + Tr [a z p] a z 
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where the Bloch vector of components (Tr [a x p] , Tr [u y p] , Tr [cr z /5]) remains on 
the unit sphere § 2 . 

The resonance assumption for the control field allows us to reduce the system 
by averaging and removing highly oscillating terms of frequency w. Consider 
the following time-dependent change of variables 



Since [a z , P] — 0, we have: 



—£ = — iAO cos(cji) 
dV ' 



—£ = — iAO cos(u)t) 
dV v ' 



o x e 



e 1 ^ 2 * a x e 



r 

-Tr 

4 



Ozt, + &z - 2Tr 



—0 =— cosM)Tr 
dt 2 K ' 



e 1 2 cr y e 1 2 £ 



Tr 



But 
and 



e 8 " 2 a *o x e l " 2 a * = e luJta *<7 x = cos(ujt)a x — sm(u)t)a y . 



- — e lu a 'o y — cos(ujt)a y + sm(Lut)a x - 

Take then the averaged system where the rapidly oscillating terms (associated 
to s'm(2uit) or cos(2w<)) are removed: 



r d AO 

dt^=- l - M 



a r X + £<r» - 2Tr 



tj (13) 



d - 7A 



Tr 



Consider now the new variables, 
C = e* 2 £e 

Then JH reads: 



+ T Tr 
4 



e 2 cr z e 



(e l 2 a z 



"(C-C) 



dt 4 



■ A0tcT x _ .Aeto 

1 2 cr„e 1 2 



X 



Tr 



e 2 CT z e 



- 2Tr 

"(C-C) 



e* 2 cr z e 



(14) 
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But we have 



■ A9t<r x ■ AOta 

Tyt 



e 2 a y e 2 = 

- A8tcr x ■ AOtcr x 



cos(A9t)ay - s\n{A9t)a z 
sm(A9t)a v + cos(A9t)a z 



Let us consider now the secular terms in (|14|) . i.e., terms without the rapidly 
oscillating factor sin(2 A9t) or cos(2 A9t). The secular term in 



Tr 



■ AOta 



■ AOta 



e 2 a z e 2 

■ AeteT x ■ AOta 



-(C-C) 



le" 2 a z e ' 2 C + (e l 2 <J z e~ 
is sum of two quantities 



2Tr 



e 2 a z e 



%(C-C) 



2Tr 
1 



-Tr 



The secular term for Tr 
ifTr 



2 c 



- 2Tr 



o- 2 C 



Tr 



e* 2 a z e 



"(C-C) 



IS 





Tr 


mc-cY 


- Tr 




Tr 


'^(C-C)' 



Thus the averaged system associated to (fT3|) reads: 



. , A(fl-fl) r . 
C = -i — — Wx,C\ 



o-.(C-C) 



<7*C + Co-» -2Tr 



o-yC Tr ct z (C-C) 



-Tr 



2Tr 


CyC 


K) 


2Tr 


CTzC 





<7 Z C 


Tr 


^(C-0 



(15) 



Consider now the following Lyapounov function: 



V((,9) = -Tr 



-Tr 



One has 
4 d 



V = Tr 



(tt^cc-C) 



Tr 



CTyC 



Tr 



^(C-C) 



Tr 



a z ( 



Tr a,(C-C) -Tr cr z (( - {) 



9 



We have (Cauchy-Schwartz inequality) 



4_d 
Ydt 



V < Tr 



Since Tr 



CyC 



+ Tr 



Tr 



Tr 



^(C-C) 



Tr 



^z(C-C) 



either 



Tr 



< 1, f t V < 0. When ^y = 0we have 



*„(C-0 



Tr 



= 0. 



This means that Tr 



c = c±- 



^(C-C) 

a x Q\ =±Tr[* x C]: 

I ± Tr [<r x C] ^ + Tr [a y (] a y + Tr [a z C] o z 



We can use here Lasalle invariance principle since £ evolves on a compact 
manifold and V is infinite when is infinite. Since £ is constant, we 
have = and thus (0 - 0)[a x ,C] = 0. Since p(0) ^ (1 ± <r x )/2, 
C(*) = C(0) ~ p(0) 7^ (1 ± a x )/2 and thus [<r x ,C] ^ 0. This implies that 
= 0. A simple inspection shows that (£_, 0) is an unstable equilibria for 
(C,0) and (C + = C,0) is an asymptotically and exponentially stable one 
(use on the first order approximation the same Lyapounov function and 
the invariance principle that proves the asymptotic stability of the first 
variation, and thus its exponentially stability). 

or (equality case for the Cauchy-Schwartz inequality) 



Tr 



<j y C =aTr[* y (], Tr 



aTr [a z (\ , 



for some real a and 



Tr 



Thus we have 



CT yC 



Tr[<r B C] 



Tr 



1. 



2 W V 



Tr[^C] 



c = c± = - 

Moreover using LaSalle invariance principle and developing dQ/dt = 0, we 
obtain = 0. One can easily see that (( — (±, = 0) are also equilibrium 
points of the system. 

However, note that we are looking for a local result. For initial state 
(Co,0o) near enough to the initial state (C)0)> we have 



V(Co,0o)<e«l. 
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As the Lyapunov function V keeps decreasing and by choosing e to be 
small enough (e < \ [\ + ^Tr [a y Q] 2 + Tr [cr 2 C] 2 ^ ), the state ((J) can 
not reach the equilibrium states (£- , 9) where we have: 



V(C-,e) = \ (l + ^Tr[a,C] 2 +Tr[a z C] 2 ) 



2 

> e. 



Concerning the equilibria (£_|_,0), we still have two situations: 

1. either Tr [cr x C] = 0, in which case = ( and we have the convergence 
as we wanted to prove in the theorem. 



2. or Tr^C] ^ which implies that yJTr [a y Q] 2 + Tr [a z C] 2 < 1. In 
this case we choose 

F(Co,#o)<e«l, 

with < e < \ (\ - -^Ti- [dyCf + Tr [cr z C] 2 ^) • As the Lyapunov 

function V keeps decreasing, the state ((, 9) can not reach the equi- 
librium states (C+, 9) where we have: 

} I 1 - » .'Tr \rr... ~ -- TV \n ( '\~ * 



V((+,6) = ^ I 1 - JTr KC] 2 + Tr [a z C] 2 ) > e. 



This asymptotic analysis based on the above Lyapounov function and Lasalle 
invariance principle shows that the steady-state (£, 9) = (C 9) of the average 
system (IT51) is locally exponentially stable. 

The existence of e and 77 results from a classical lemma concerning the aver- 
aging techniques (cf. [9], Page 333, Theorem 8.3): if the average system admits 
an exponentially stable steady-state that is also an equilibrium of the origi- 
nal system, then this steady-state is also exponentially stable for the original 
system. 

Consider first (fT5l): 



• It admits (£, 9) as exponentially stable steady-state. 

• it corresponds to the averaging of (fT*!)) when A9 ^> T, 7 A 

Thus the state ((,9) of (|14[) converge exponentially towards (£, 9). Since 

£ — £ and 9 — 9 converge exponentially towards 0. when to 3> AO, a similar 
argument between (fTTj) and (fl3l) yields to the local exponential convergence of 
p — p and — towards 0. 

□ 
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Remark 3. Theorem [7] provides a local convergence result for the estimator 
equations given by (jlip . The simulations, however, show a much stronger global 
convergence behavior. Notice that we cannot prove directly that (fT3"|) is stable, 
although the following function 



is decreasing in average since 



■Tr 



TV 
dt 



2^ 



I - Tr 



2 H 

1 



AOTr 



The first term of the left hand side is always negative whereas the second one 
admits a zero average. Thus in average, V is a decreasing function of time. 



3 The general case 

In this section, we extend the above identification algorithm to the general 
case of a multi-level system. Before going through this extension, we need 
to address the non-trivial identifiability problem for the multi-level case. The 
AppcndixlAl (based on the result of |15j ) provides sufficient assumptions ensuring 
the identifiability. 

3.1 Formal extension 

Consider the iV-levels system, (\j))i<j<N, described by the density matrix p 
that obey the following dynamics (we assume here the assumptions A1,A2 and 
A3 of the Appendix lAl to be satisfied): 

[j t P = -i[H + u{t)ii,p] (16) 

where 

• H = ^2f =1 tjJj \j) (j\ is the free Hamiltonian with luj real and satisfying 
\u>i — u)k\ 7^ \<^v ~ for an y distinct couples (I, k) and (/', k'). 

• p, = J2i<i<k<N ®ik (\k) (l\ + \l) (k\) where 6^ are the parameters to iden- 
tify 

• the electromagnetic field is represented by the scalar input u(i) eM 
We assume that 

Vj (t) = Tr[P jP (t)], P 3 = \j)(]\, j = l,2,...,N 

are the measured outputs. The goal is to estimate the coefficient The u>j 
arc known. 
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The estimator (fTTj) admits then the following generalization (1 < I < k < N): 



N 



-p = -i[H + u{t)p, p] (%(*) - Tr [Pjp]) {P 3 p + pP s - 2Tr [Pjp] p) 



3=1 

N 



±9 lk = -i llk u(t) Tr [P 3 (k\ + \k) (l\ , p}} ( Vj (t) - Tr [P jP ]) 

3=1 

(17) 

where p = J2i<i<k<N^k (10 ( fc l + (^1)' r > and Hk > are design pa- 
rameters. Notice that p remains a projector if its initial condition is also a 
projector. 

Set u(t) = J2i<i<k<N Aik cos(ujikt) where ui k = ui~uj k and A lk is a constant 
amplitude. Let us compute, at least formally, the averaged system associated 
to dTTJ) when 

< T < Ai k 9i k < uji k , yik<.0i k . 
Consider the "Pauli matrices" associated to the transition between I and k: 

<7 l x k = \l)(k\ + \k)(l\, a l y k = -i\l){k\+i\k){l\ 

a? = P l -P k = \l) (l\ - \k) (k\ , I lk =P l+ P k = \l) (l\ + \k) (k\ 

For each / ^ k, we have the usual relations: 

(~lk\2 jlk J.kJ,k -Ik 

{<j x ) = 1 , <r x a y = i<j z , ... 

jjk I jfc 

For each j and k ^ j, Pj = I n the sequel, we use the shortcut notation 

^ ;fc that stands for J2i<i<k<N- Thus we have 

u = A ifc cos(wifci), /i = 0jkoi fc . 

Notice that when j ^ I and j ^ k, Tr [P 3 [a lk , p}] = 0. When j = I, set 
P = (P k +<J l z k )/2 to find 

Tr[P l [a l x k ,p]]=iTr[ ( r l y k p] 
When j = k we have similarly 

Tr[P fc [^ fc ,/5]] = -tU: . 

Thus (fT7| reads (remember that Pi — P k = & lk )- 
d 



, P = ~i[H,p] -iJ2J2 Al ' k '® lk cos{u vk d)[a x k ,p] 

Ik l'k' 

N 

+ r £ Tr ^ (p - p)] (P,p + pP j - 2Tr [P;p] p) 



Ik 



dt 



3=1 

7ife (E <*»K**) J Tr [a lk p] (Tr [af (p - p)] ) 



sl'k' 
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In the inter-action frame, £ = e pe and £ = e pe , we have 



; 2J M'k'Olk cos(u>i'k't) 



AHt JLk „ — iHt 

e a„ e , 



;fc /'ft' 
N 



+ r^Tr[p j (£-0 



j'ft' 



iff* (fc -iHtt 

e cr e £ 



Tr 



since [Pj,H] = 0. Simple computations show 



iHt Ik -iHt _ iuji k ta lk Ik 



e a, r e 



a 1 * = coB(w lk t)o£ - sin(wifct)a^ 



jft 



and 



e iHt a lk e~ iHt = e^IV 



„ - sin(w zfc i)CT^ fc + cos(cJikt)a l y K 

For (/,£;) ^ (l',k'), \u>ik\ 7^ l^i'fc'l- Thus resonant terms come only from (Z,fc) 
(l',k'). The "rotating wave approximation" of (jTTJ) reads: 



/ft 

N 



A 



lk u lk 



r^Tr[F,(e-0 



P^ + ^-2Tr 



(18) 



6a= TU^ Tc [o**i]'I t [o*it-£) 



Notice that, instead of using (fl~7|) as estimator, one can use in practice such 
averaged filter where the large transition frequencies u>ik are removed: 



dt 



EAikOik 
2 

Ik 
N 

Ilk Alk , 



x )? 



Tr 



AC 



P,£ + C-Pj - 2Tr 



AC 



(19) 



- Tr 



Let us now generalize the heuristic argument of remark[3]for the stability of (|13p . 
Consider the following function 



1 W 2 



2(^;fc — Oik) 2 
Ilk 



(20) 
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One can easily see that 



dV_ 
~dt 



EAlkOlk ^ 
2 1 

Ik 



- 2r EE Tr [ p ^ 



Tr 



(21) 



While the second term in (j2Tj) is obviously negative, the first term has no reason 
to be negative. However, we will show by a formal argument that (considering 
some appropriate assumption concerning the Rabi frequencies) this term can be 
averaged to zero and thus can be neglected. 

In this aim, consider the real effective Hamitonian: 



H, 



eff 



E 

Ik 



AikOik i k 

r 0% , 



and diagonalize it as follows: 

H e f f =EflClE, Sl = diag(n 1 ,...,n N ), E lk € 



VI, k. 



where {fij}jLi are Rabi frequencies of the system. From now on, we will assume 
that these Rabi frequencies are non-degenerate (fi m ^ fl n for m ^ n) and 
moreover that T <C Aq and 7; A; <C An, where Aq 

Now, in analogy with the 2-level case, consider the unitary transformation 

C=U*EZEftU, C = U^EiE^U, 

where U(t) = exp(— it£l). Under such a transformation £ is trivially constant 
( = E^qE^ . Furthermore, this transformation also removes the highly oscillating 
part of £, (\§ik -Oik\ < 6 m n and T < Ai k 8i k for all I, k, m, n): 

= _.^ A lk (O lk -O lk ) [rfEjhtfjj^ 

Ik 

N 

+T^Tr ^EPjE^UiC - O] (u^EPjE^U^ + (U*'EP 3 E*<U - 2Tr [t/ t £P ) .B t l7c] C 

j=i 

Now let us develop the terms in the first part of (|2"Tj) using this unitary trans- 
formation: 

Tr 
Tr 



[u f Ecrl k E f U(( - O] Tr \u ] Ea l y k E ] U{C - C)] 



J2 eX P (i(fi r - Sl s )t)(E H E sl - E rk E sk )(Csr ~ £sr) + ^(^r! - E 2 rk ){C,rr ~ Crr) 



^exp(i(fi r - Q s )t)(E rl E sk - E rk E 3l )(( ar - ( sr ) 
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Developing and removing the highly oscillating terms of frequencies An, we find 

2_.(E r lE s l — E r kE s k)(E r kEsl — E r lE s k)\^sr — Csr| 2 = 

- ^>~]{(E r iE s i — E rk E sk )(E r iE 3k — E r kE„i)\Crs — Crs\ 2 + 
^ _, 

(E r lE 3 l — ErkEsk)(E a lE r k — E a kE r l)\Csr ~ Or) 2 \ = 0, 

where we have broken the sum into two parts by symmetrizing with respect to 
the indices r and s. 

Even though this argument does not prove the convergence of the estimator 
for the multi-level system, it gives a strong reason for it to be efficient. The 
simulations of the next section show that this is effectively the case. 

3.2 Simulations 

Let us check out the performance of this algorithm on two other test cases: first 
on a 3- level system and next on the 4- level system also considered in [TS] . 
The first test case is given by 



/0 0\ 
H= 1 
\0 3/ 







1.3 


and 


[i = 1 1.3 









-1.5 



We consider a laser field resonant with all the transition frequencies of the 
system: 

u(t) = A (sin(t) + sin(2t) + sin(3t)) , 

where the laser amplitude A here is chosen to be 0.1. This allows us to use the 
rotating wave approximation (averaging) and eliminate the Hamiltonian Hq in 
order to obtain an effective Hamiltonian H e ff. The Figures [3] and |4] illustrate the 
simulations of the filter equation (fl~7|) where the following constants have been 
used: 7 = 1, V = .05, A = 0.1. The simulation time T — 1600 represents around 
250 times the largest period of the system's Hamiltonian H . This, however, 
represents only about 25 times the largest period of the effective Hamiltonian 
H e ff. Furthermore, the parameter estimator /i, the state estimator p = 
and the real state p = "J"!^ are initialized at: 

/ 1.2 .9 \ A/Vi4\ A/\/30\ 

A(0)= 1.2 -1.7 , *(0)= 2/VTi , *(0)= 2/V30 • 
V-9 -1.7 / \3/VuJ \5/V30j 

As one can easily see, the estimator /t ends up giving a really good approximation 
of the true dipole moment in a completely reasonable time. 

Let us now consider the 4-level test case also considered in [TS] . This permits 
us to have a comparison between the algorithm provided in this paper and the 
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Figure 3: This figure illustrates the evolution of the estimator /i; the estimator 
gives a very good approximation of the real dipole moment /x. 
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Figure 4: The positive definite function V defined in (|20| with 7^ = .5; the 
function decreases in average and ends up converging toward 0. 



numerical one presented in 



H = 



I 0.0833 
-0.0038 
-0.0087 

V 0.0041 



-0.0038 
0.0647 
0.0083 
0.0038 



15] . The physical system is given by: 

-0.0087 0.0041 \ 

0.0083 0.0038 

0.0036 -0.0076 

-0.0076 0.0357 / 



f° 


5 


-1 





5 





6 


-1.5 


-1 


6 





7 


\o 


-1.5 


7 






Diagonalizing the matrix H yields 

/0 




H = P D P- 



where 



D = 









0.0365 








0.0651 








0.0857/ 



V 



/o 1 

-1 

1 -1 



-1 

1 1 

-1 

1 0/ 



P = exp(n 
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is an anti-Hcrmitian matrix. We consider a laser field resonant with all the 
transition frequencies of the system: 

4 

u(t)=A 5^53sin((A,-A fc )t), 

1=1 k<l 

where Xj represents the j'th eigenvalue of H and the laser amplitude A is chosen 
to be 0.01. 



4. Si 
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Figure 5: This figure illustrates the evolution of the estimator (1 towards; the 
estimator gives a very good approximation of the real dipole moment /x. 




Figure 6: The positive definite function V defined in (|20| with 7^ = .5; the 
function decreases in average and ends up converging toward 0. 

The Figures [5] and [5] illustrate the simulations of the filter equation (|17[) 
where the following constants have been used: 7 = 0.5, T = 1, A = 0.01. The 
simulation time T = le + 05 represents around 320 times the largest period of 
the system's Hamiltonian H . This represents about 650 times the largest period 
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of the effective Hamiltonian H e ff. Furthermore, the parameter estimator /t, the 
state estimator p — and the real state p = are initialized at: 



A(o) = 



/ 6 -1.5 .05^ 

6 7 -2 

-1.5 7 6 

\ .05 -2 6 / 



tf(0) = 



fi/vm 

2/V30 
3/V30 

WvW 



*(0) = 



N 2 \ 

1/2 
1/2 

v/y 



As one can easily see, the estimator /t ends up giving a really good approximation 
of the true dipole moment. 

4 Robustness 




Figure 7: Here we consider the noised measurement of the form (|2"2"|) ; the 
first figure illustrates the robustness, with respect to the uncertainties, in the 
evolution of the parameter estimator 8; the second one shows the robustness in 
the evolution of the error term e(t). 
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Figure 8: The robustness in the evolution of the Lyapunov type function V with 
7 =.l. 



The laboratory noises are always present and are not negligible. These noises 
affect both the output result and the laser field. Moreover the delay in reading 
the laboratory output results are essential and must be taken into account in a 
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Figure 9: Here we consider the noised control field of (|22[l ; the first figure 
illustrates the robustness, with respect to the uncertainties, in the evolution of 
the parameter estimator 9; the second one shows the robustness in the evolution 
of the error term e(t). 
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Figure 10: The robustness in the evolution of the Lyapunov type function V 
with 7 = .1. 



faithful model. In this section, we study the robustness of the algorithm with re- 
spect to all these uncertainties through a number of simulations. Concerning the 
measurement output results, three kind of uncertainties can be admitted: a de- 
lay in reading the output result, a small additional constant gain and additional 
non-correlated noises. The simulations show that the identification algorithm, 
presented above, is robust with respect to all these uncertainties. The simu- 
lations of Figures [7] and [5] show this fact for the 2-level system of Section [51 
Here we have considered a delay of 0.3 (about 1/2 of the natural period of the 
system) in reading the measurement results. Moreover, we have added small 
additional constant gains and additional non-correlated gaussian noises: 

y{t) = Tr [P p(t - 0.3)] + 0.06 + 0.07 w, 

where w has a standard normal distribution. Other simulation parameters are 
fixed exactly as in the Section [2] 

Regarding the control input, we consider two kind of uncertainties: a small 
additive constant gain for the amplitude A of the laser field and an additional 
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gaussian noise for the laser field. The simulations of Figures |H] and [TUl on the 
2-level system of Section [2j show the robustness of the algorithm with respect 
to these uncertainties. Here, we have assumed that the laboratory laser field is 
noised as follows: 

u(t) = (A + 0.03) sin(t) + 0.07 w, (22) 

where A — 1 as in Section[2]and w is a normal distribution. Similar simulations 
concerning the systems of higher dimensions represent the same kind of robust 
behavior. 

Remark 4. One might consider additional uncertainties concerning the fre- 
quencies of the laser field (e.g. small additive constant gains in the laser fre- 
quencies). As it has been discussed in Remark^ a more realistic estimator in 
the settings of our paper is given by the first averaged filter (|13p ( (|19p in the 
case of a multi-level system). The Bohr frequencies of the system do not appear 
in this estimator. Therefore, one can easily check that this averaged estimator 
represents a robust behavior with respect to the uncertainties in the laser fre- 
quencies. One only needs these frequencies to be near enough to the transition 
frequencies of the system (uj r — u <C A9 in the case of the 2-level system). 

5 Conclusion 

In this paper, we propose an observer-based method for the Hamiltonian iden- 
tification of a quantum system. An intuitive observer ^ has been considered 
for the Schrodinger equation (JTJ) and has been developed and extended to give 
an estimate of the unknown parameters of the system. The convergence of this 
method is completely analyzed for a 2-level case. The multi-level cases have 
been addressed using heuristic arguments. Various simulations in different di- 
mensions illustrate the relevance of the technique for these multi-dimensional 
systems. Finally the robustness of the design with respect to different uncer- 
tainties and noises is addressed by simulations on the 2-level case. Similar 
robustness results can be noted for multi-dimensional systems. 

In RemarkUl it has been noted that replacing the estimator (|17p with the first 
averaged version (fTl?)) , increases considerably the robustness of the identification 
result with respect to the frequency uncertainties. 

Such averaged filter represent even more advantages whenever the settings 
considered in the paper are valid. In particular, one increases considerably the 
robustness with respect to the delay in the measurement. Indeed, this delay only 
needs to be much smaller than the shortest Rabi period of the system. Secondly, 
the non-degeneracy assumption for the Rabi transitions An may be removed 
using a slow modulation of the amplitudes A^. Finally, one does not really 
need to have access to the continuous measurement results yj(t) (which is lots 
of information to be asked in the laboratory settings). In fact, one only needs 
samples on the output signal with frequencies much higher than the larger Rabi 
frequency. All these advantages seem to highly privilege the use of the averaged 
estimator (|T5)) . 
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A Identifiability 

In this appendix, we present the mathematical framework in which the identifi- 
cation problem can be considered. Moreover, we review briefly the former work 
on the identifiability of the considered system. A well-posedness result which 
allows us to consider the identification problem in Section [3] will be announced. 
The goal is to identify H = Hq + V or/and fi in system 

i^9 = (H + V + u{t)n)% *| i=0 = *o, ||*o||w = l, (23) 

when laboratory measurements on some physical observables are provided. In [15j , 
two different settings have been considered in order to characterize the identifi- 
ability of such a system: 

(51) The Hamiltonian H is known and the goal is to identify the dipole mo- 
ment fi. The so-called populations along the eigenstates <pi, i.e. pi = 
| (<pi, \ 2 ,i = l,2,...,N are measured for all instants t > 0. This is 
performed with as many control amplitudes u(t) as required. 

(52) Neither the potential V nor the dipole moment fi are known and the goal 
is to identify them. Note that, by identifying H we mean identifying V, 
as Hq is readily known. The eigenvalues of the Hamiltonian H = Hq + V 
are also assumed to be known (this assumption is relevant in practice, 
see Remark [5]). Here we measure the populations pi along the states of a 
canonical basis {ei}fL 1 : Pi = \ (e,-, V?(t)) \ 2 ,i — 1,2,..., A for all instants 
t > and all control amplitudes u(t). 

Remark 5. It is relevant in practice to assume that the eigenvalues of the 
internal Hamiltonian H = Hq + V are known. In fact the classical spectroscopy 
allows for identifying the eigenvalues of the Hamiltonian and discriminating 
between two systems that do not share the same ones. In fact spectroscopy only 
gives eigenvalue differences (transition frequencies), not the absolute values. The 
overall unknown additive factor is not seen by the measurements and has no 
impact on the identification result. 

In this paper, we have only considered the first setting. An extension of the 
technique to the second setting remains to be done in future work. However, 
[15] provides an identifiability result for this second setting as well. 

Here we announce the identifiability result of |15j concerning the first setting. 
For a result in the second setting and also the proof of the result for the first 
setting, we refer to [15] . 



24 



Theorem 2. Suppose that there exist two dipole moments [i\ and fj,2, giving 
rise to two evolving states and ^2 respectively solving 

i*i = + (24) 
ii> 2 = (H + u(t)M 2 )*2, (25) 

that produce identical observations for all t > and all fields u(t): 

|<*i(f),<&) | 2 = |(* 2 (*),^)| 2 t= 1,2, ...,7V. (26) 

TTien under assumptions 

(Al) Equation (|24[) is wavefunction controllable \25jj : 

(A2) TVie transitions of the Hamiltonian H are non-degenerate: A^ — 7^ 
A 42 - A J2 /or (*i,ji) ^ (i 2 , 32) WB; 

(A3) TTie diagonal part of the dipole moments [i\ and 112, when written in 
the eigenbasis of the Hamiltonian H, is zero: {4>\i^i \4>)i = (^IjA^ \4>)i — 
0,*= 1,2,.. .,7V; 

i/ie too dipole moments are equal within some phase factors {cti}f =l C R such 
that: 

Vi,j = l,2,...,N, =e i ^-^{n 2 ) ij . (27) 

Fore more details and remarks concerning the assumptions and the result of 
this theorem we refer to [15] . 
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